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Abstract. In this work we develop and analyze a mathematical model of bi¬ 
ological control to prevent or attenuate the explosive increase of an invasive 
species population in a three-species food chain. We allow for finite time blow¬ 
up in the model as a mathematical construct to mimic the explosive increase in 
population, enabling the species to reach “disastrous” levels, in a finite time. 
We next propose various controls to drive down the invasive population growth 
and, in certain cases, eliminate blow-up. The controls avoid chemical treat¬ 
ments and/or natural enemy introduction, thus eliminating various non-target 
effects associated with such classical methods. We refer to these new con¬ 
trols as “ecological damping”, as their inclusion dampens the invasive species 
population growth. Further, we improve prior results on the regularity and 
Turing instability of the three-species model that were derived in [43]. Lastly, 
we confirm the existence of spatio-temporal chaos. 


1. Introduction 

Population dynamics are a fundamental aspect of many biological processes. In 
this paper, we introduce and investigate a mathematical model for the popula¬ 
tion dynamics of an invasive species in a three-species food chain model. Exotic 
species are defined as any species, capable of propagating themselves into a non¬ 
native environment. If the propagating species is able to establish a self sustained 
population in this new environment, it is formally defined as invasive. The survival 
and competitiveness of a species depends intrinsically on an individual’s fitness and 
ability to assimilate limited resources. Often invasive species possess the ability to 
dominate a particular resource. This allows them to expand their range via out- 
competing other native species. In the United States damages caused by invasive 
species to agriculture, forests, fisheries and businesses, have been estimated to be 
$120 billion a year [47]. In the words of Daniel Simberloff: “ Invasive species are a 
greater threat to native biodiversity than pollution, harvest, and disease combined .” 
[54] Therefore understanding and subsequently attenuating the spread of invasive 
species is an important and practical problem in spatial ecology and much work has 
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been devoted to this issue [1, 4, 10, 38, 40, 53]. More recently, the spread of natural 
and invasive species by nonrandom dispersal, say due to competitive pressures, is of 
great interest [2, 35]. However, there has been less focus, on the actual eradication 
of an invasive species, once it has already invaded. This is perhaps a harder prob¬ 
lem. In the words of Mark Lewis: “Once the invasive species are well established, 
there is not a lot you can do to control them”. [66] It is needless to say however, 
that in many ecosystems around the world, invasion has already taken place! Some 
prominent examples in the US, are the invasion by the Burmese python in south¬ 
ern regions of the United States, with climatic factors supporting their spread to a 
third of the United States [ ]. The sea lamprey and round goby have invaded the 
Great Lakes region in the northern United States and Canada [7]. These species 
have caused a severe decline in lake trout and other indigenous fish populations. 
Lastly the Zebra Mussel has invaded many US and Canadian waterways causing 
large scale losses to the hydropower industry [39]. 

Another factor attributed to the increase of an invasive population, is that the 
environment may turn favorable for the invasive species in question, while becoming 
unfavorable for its competitors or natural enemies. In such situations, the popula¬ 
tion of the invasive species may rapidly increase. This is defined as an outbreak in 
population dynamics [6] . These rapid changes tend to destabilize an ecosystem and 
pose a threat to the natural environment. As an illustration, in the European Alps 
certain environmental conditions have enabled the population of the larch budmoth 
to become large enough that entire forests have become defoliated [36]. In most 
ecological landscapes, due to exogenous factors, one always encounters an invasive 
species and an invaded species. If the density of the former, be it invasive, disease 
causing, an agricultural pest, defoliator or other, undergoes a rapid transition to a 
high level in population, the results can be catastrophic both for local and nonlocal 
populations. 

Biological and chemical controls are an adopted strategy to limit invasive pop¬ 
ulations [65]. Chemical controls are most often based on direct methods, via the 
use of pesticides [65]. Biological control comprises of essentially releasing natural 
enemies of the invasive species/pest against it. These can be in the form of preda¬ 
tors, parasitoids, pathogens or combinations thereof [65]. There are many problems 
with these approaches. For example, a local eradication effort was made by USGS 
through a mass scale poisoning of fish in order to prevent the asian carp (an in¬ 
vasive fish species) from entering the Chicago Sanitary and Ship Canal [30]. The 
hope was to protect the fishing interests of the region. However, among the tens of 
thousands of dead fish, biologists found only one asian carp. Thus chemical control 
is not an exhaustive strategy. However, biological controls are also not without 
its share of problems. In fact, sometimes the introduced species might attack a 
variety of species, other than those it was released to control. This phenomena is 
referred to as a non-target effect [16, 34], and is common in natural enemies with a 
broad host range. For example, the cane toad was introduced in Australia in 1935 
to control the cane beetle. However, the toad seemingly attacked everything else 
but its primary target [46]! In addition, the toad is highly poisonous and there¬ 
fore predators shy away from eating it. This has enabled the toad population to 
grow virtually unchecked and is today considered one of Australia’s worst invasive 
species [55]. In studies of biological control in the United States estimate that when 
parasitoids are released as biological controls that 16 % of the introduced species 
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will attack non-targets, in Canada these numbers are estimated as high as 37.5% 
[16]. In practise it is quite difficult to accurately predict these numbers. The cur¬ 
rent drawbacks make it clear that alternative controls are necessary. Furthermore, 
that modeling of alternative controls is important to validate the effectiveness of 
a management strategy that hopes to avoid non-target effects. Such modeling is 
essential to access and predict the best controls to employ, so that the harmful pop¬ 
ulation will decrease to low and manageable levels. This then gives us confidence 
to devise actual field trials. We should note that in practice actual eradication is 
rarely achieved. 

Thus there are clear questions that motivate this research: 

(1) How does one define a “high” level for a population, and further, how well 
does an introduced control actually work, at various high levels? 

(2) Is it possible to design controls that avoid chemicals/pesticides/natural 
enemy introduction, and are still successful? 

This paper addresses these questions through the investigation of a mathematical 
model that: 

(1) Blows-up in finite time. Given a mathematical model for a nonlinear pro¬ 
cess, say through a partial differential equation (PDE), one says finite time 
blow up occurs if 

lim ||r||x -> oo, 

t —> T * <oo 

where X is a certain function space with a norm || • ||, r is the solution to the 
PDE in question, and T* is the blow up time. Therefore “highest” level is 
equated with blow up, and the population passes through every conceivable 
high level of population as it approaches infinity. 

(2) Incorporates certain controls that avoid chemicals/pesticides/natural en¬ 
emy introduction. The controls we examine are: 

(a) The primary food source of the invasive species is protected through 
spatial refuges. The regions that offer protection are called prey refuges 
and may be the result of human intervention or natural byproducts, 
such as improved camouflage. 

(b) An overcrowding term is introduced to model the movement or disper¬ 
sion from high concentrations of the invasive species. Densely popu¬ 
lated regions have increased intraspecific competition and interference 
which cause an increase in the dispersal of the invasive species. This 
is an improvement to current mathematical models and will be seen 
to be beneficial if used a control. 

(c) We introduce role reversing mechanisms, where the role of the primary 
food source of the invasive species, and the prey of this food source, 
is switched in the open area (the area without refuge). This models 
situations where the topography provides competitive advantages to 
certain species. It will be seen that this also is beneficial if used a 
control. In effect, this uses the current ecosystem and by modifying 
the landscape a natural predator in the environment has an advantage 
in key areas. Hence, the invasive species population will adversely be 
effected. It can also be thought of as introducing a competitor of the 
invasive species, to compete with it for its prey. 
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Remark 1 . Note, none of the above rely on enemy release to predate on the 
invasive species, or a parasite or pathogen release to infect the invasive species. 
Thus potential non-target effects due to such release can be avoided. 

In the literature finite time blow up is also referred to as an explosive instability 
[58] . There is a rich history of blow up problems in PDE theory and its interpre¬ 
tations in physical phenomenon. For example, this feature is seen in models of 
thermal runaway, fracture and shock formation, and combustion processes. Thus 
blow up may be interpreted as the failure of certain constitutive materials leading to 
gradient catastrophe or fracture, it may be interpreted as an uncontrolled feedback 
loop such as in thermal runaway, leading to explosion. It might also be interpreted 
as a sudden change in physical quantities such as pressure or temperature such 
as during a shock or in the ignition process. The interested reader is referred to 
[48, 58]. Blow up in population dynamics is usually interpreted as excessively high 
concentrations in small regions of space, such as seen in chemotaxis problems [26]. 
Our goal in the current manuscript is to bring yet another interpretation of blow 
up to population dynamics, that is one where we equate such an excessive con¬ 
centration or “blow up” of an invasive population with disaster for the ecosystem. 
Furthermore, it is to devise controls that avoid non-target effects and yet reduce 
the invasive population, before the critical blow up time. 

In the following, the norms in the spaces L p (fl), L°°(n) and C (fl) are respec¬ 
tively denoted by 



n 


In addition, the constants C, C\ and C 2 may change between subsequent lines, as 
the analysis permits, and even in the same line if so required. 

The current manuscript is organised as follows. In section 2 we formulate the 
spatially explicit model that we consider. In section 3 we describe in detail the mod¬ 
eling of the control mechanisms that we propose, and term “ecological damping”. 
Here we make three key conjectures 1, 2 and 3, concerning our control mechanisms. 
Section 4 is devoted to some analytical results given via lemma 4.3, theorem 4.4 and 
4.5. Section 5 is where we explain our numerical approximations and test conjec¬ 
tures 1, 2 and 3 numerically. In section 6 we investigate spatio-temporal dynamics 
in the model. We investigate the effect of overcrowding on Turing patterns, and we 
also confirm spatio-temporal chaos in the model. Lastly we offer some concluding 
remarks and discuss future directions in section 7. 


2. Model Formulation 


A three species food chain model is developed, where the top predator, denoted 
as r, is invasive. In our model, r may blow up in finite time. Although populations 
cannot reach infinite values in finite time, they can grow rapidly [6]. The blow 
up event indicates that the invasive population has reached “an extremely high” 
and uncontrollable level. Naturally, that level occurs prior to the blow up time. 
Therefore, as ||r|| approaches infinity in finite time, it passes through every con¬ 
ceivable “high” level. The blow up time, T*, is viewed as the “disaster” time. We 
investigate mechanisms that attempt to lower and control the targeted population 
before time T*. This modeling approach has distinct advantages: 
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(1) There is no ambiguity as to what is a disastrous high level of population. 

(2) There is a clear demarcation between when or if the disaster occurs. 

(3) The controls that are proposed do not rely on a direct attack on the invasive 
species, as is the traditional approach, rather we attempt to control the food 
source of r. This will avoid possible nontarget effects. 

(4) The model provides a useful predictive tool that can be tuned and estab¬ 
lished through data, in various ecological settings. 

(5) Mathematical models are advantageous in many situations due to their cost- 
effectiveness and versatility. Of course obtaining an analytical solution for 
a nonlinear model is virtually impossible, outside of special cases. However, 
a numerical approximation can be developed to accurately investigate the 
role and effect our controls have on the blow up behavior. 

Suppose an invasive species r has invaded a certain habitat and it has become 
the top predator in a three species food chain. Hence, r predates on a middle 
predator v, which in turn predates on a prey u. A temporal model is given for the 
species interaction, namely 


( 1 ) 

( 2 ) 

( 3 ) 


dr 

dt 

dv 

dt 

du 

dt 


= cr — w 3 


v + D 3 ’ 


— = -a 2 v + wi 


u + D i 


W 2 


V + D 2 


VAjUj Q 

— = a\u — b 2 u — w o 


u + Dq 


The spatial dependence is included via diffusion, 

r 2 

d t r = d 3 Ar + cr 2 - ^ = h(u,v,r), 


( 4 ) 

( 5 ) 

( 6 ) 


v + D 3 

_ , . uv 

o t v = d 2 Av — a 2 v + w\ -— 

u + Di 

dtu = d\ Au + a\U — b 2 u 2 — Wg 


— w 2 


vr 


v + D 2 


= g(u,v,r), 


u + D( 


= f(u,v,r), 


defined on M + x fh Here H C N = 1,2 and A is the one or two dimensional 
Laplacian operator. We define x to be the spatial coordinate vector in one or two 
dimensions. The parameters d\. d 2 and d 3 are positive diffusion coefficients. Neu¬ 
mann boundary conditions are specified on the boundary. The initial populations 
are given as 


u(0,x) = u 0 (x), u(0, x) = v 0 (x), r(0, x) = r 0 (x) in fi, 


are assumed to be nonnegative and uniformly bounded on Q. 

There are various parameters in the model: a±, a 2 , b 2 , wq, w i, w 2 , «; 3 , c, Dq, 
D i, D 2 , and D 3 are all positive constants. Their definitions are as follows: a 3 is 
the growth rate of prey u; a 2 measures the rate at which v dies out when there is 
no u to prey on and no r; Wi is the maximum value that the per-capita rate can 
attain; Dq and D\ measure the level of protection provided by the environment 
to the prey; b 2 is a measure of the competition among prey, u; D 2 is the value of 
v at which its per capita removal rate becomes w 2 /2 ; D 3 represents the loss in r 
due to the lack of its favorite food, v; c describes the growth rate of r via sexual 
reproduction. 
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These models offer rich dynamics and were originally proposed in [60, 61] in 
order to explain why chaos has rarely been observed in natural populations of 
three interacting species. The model stems from the Leslie-Gower formulation and 
considers interactions between a generalist top predator, specialist middle predator 
and prey. The study of these models have generated much research [18, 32, 33, 
41, 43, 62, 63, 6' ]. An interesting fact is if c > the spatially independent 
and spatially dependent models are easily seen to blow up in finite time [42]. The 
spatially dependent model offers further rich dynamics, in particular the possibility 
of Turing instabilities and non Turing patterns [32, 43]. Nevertheless, to avoid blow 
up it appears that one must restrict c < where K < 1. This was established 

for the spatially independent model in [3] and is offered here for convenience: 


Theorem 2.1. Consider the model (l)-(3). Under the assumption that 

(7] I _ w 0 b 2 P 3 _ \ ws_ 

(ai + j + w 0 b 2 D 3 J D 3 

all non-negative solutions (i.e. solutions initiating in of (l)-(3) are uniformly 
bounded forward in time and they eventually enter an attracting set A. 


Remark 2. We have recently shown the above result to be false in the ODE and 
PDE cases. That is, (4)-(6) may blow up in finite time, even under (7) provided 
the initial data is large enough [44]. It is clear that even if (7) is maintained r can 

blow up. This becomes more evident if we consider the coefficient c-"A on r 2 

and the fact that if the fecundity of r is large compared to then blow up will 

occur. This may happen in situations where: 

(1) There is an abundance of v, 

(2) r possesses certain abilities to out compete other species, and harvest 
enough v, 

(3) The environment has turned favorable for r and unfavorable for its natural 
enemies or competitors. Thus it can harvest v unchecked. 


Remark 3. If no intervention is made we can envision r growing to disastrous 
levels with adequate initial resources. The blow up time T* is viewed as the point 
of disaster in an ecosystem. Thus one is interested in controlling r via the use of 
biological controls, before the critical time T*. 

Remark 4. These observations motivate an interesting question. Assume that 
both (l)-(3) and its spatially explicit form (4)-(6) blow up in finite time for c < ^ 
for a given initial condition. Can we modify (4)-(6), via introducing certain controls, 
so that now there is no blow up, for the same initial condition? 


3. Delay and Removal of blow up via “Ecological Damping” 

It is clear that controlling the population of an invasive specie is advantageous, 
and most often necessary. However the avenues for which this is possible, whilst 
avoiding non-target effects, is not clear. Here, we propose new controls that may 
delay or even remove blow-up in the invasive population. We refer to these controls 
as “ecological damping”, akin to damping forces such as friction in physical systems, 
that add stability to a system. The crux of our idea is to use prey refuges, in 
conjunction with role reversal and overcrowding effects. There is a vast literature on 
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prey refuge, spatial refuges as well as role reversal in the literature. The interested 
reader is referred to [12, 15, 29]. However, to the best of our knowledge these have 
not been proposed as control mechanisms for invasive species. 
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FIGURE 1 . An example of a prey refuge function bi(x) = 
i (i — tanh( ^)) 

in the one dimensional case. This function decreases 
monotonically throughout its domain. The range is 0 < bi(x) <1. If 
bi(x) = 1 then the prey v is protected while if bi(x) = 0 then v is 
unprotected. 


3.1. Prey Refuge: Modified Model I. We consider modeling a prey refuge. 
Consider the continuous function 6i(x). The region where &i(x) = 1 or sufficiently 
close to one is defined as a prey refuge domain or patch. We call the region where 
&i(x) = 0 or sufficiently close to zero an open area, this is the region where r can 
predate on v. In Figure 1 we see a sharp gradient between the prey refuge domain 
and the open area. 

The inclusion of a prey refuge influences the equations for r and v , namely, 

r 2 

(8) r t = d 3 Ar + cr 2 -w 3 - 


’ (1 - bi))v + D 3 


(9) 


v t = d 2 Av — a 2 v + wi 


u + D i 


- w 2 (1 - bi) 


A 


posed on a bounded domain in one or two dimensions. Neumann boundary condi¬ 
tions are specified. The equation for u may be found in (6). 

The introduction of the prey refuge creates regions where it is impossible for 
r to predate on v. Notice, that if the entire spatial domain is considered a prey 
refuge, that is, &i(x) = 1 Vx G Cl, then the coefficient of r 2 will depend on the 
sign of c < jA. If this is negative, then the invasive species dies off. Likewise, in 
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the absence of a prey refuge, that is, &i(x) = 0 Vx G H the equations for r and 
v collapse to our previous ones, that is, (4) and (5), respectively. In such a case, 
it is know that blow up may still occur for sufficiently large enough data even if 
c < ^|. This is because the coefficient of r 2 may still be positive, as one always 
has v +b 3 < c < Zv The introduction of the refuge forces the coefficient of r 2 to 
change sign between the refuge and the open area. In the literature such problems 
are referred to as indefinite parabolic problems. The word indefinite refers to the 
sign of the coefficient being indefinite. Although there is a vast literature on such 
problems for single species models [19, 49], there is far less work for systems. There 
is also a large amount of literature for such switching mechanisms incorporated to 
understand competitive systems [20, 21], particularly in the vein of human economic 
progress. 

This indefinite parabolic problem motivates a collection of questions: If blow up 
occurs for particular parameters in the case &i(x) = 0, will a prey refuge prevent 
blow up? How does this depend on the geometry of the refuge? What is the critical 
size or shape of a refuge that prevents blow up from occurring? Does this depend 
on the size of the initial condition? In the situation that blow up still persists, how 
is the blow up time affected? Are these results influenced by multiple refuges? 

In either case, we make the conjecture: 

Conjecture 1. Consider the three species food chain model (4)-(6), a set of pa¬ 
rameters with c < and an initial condition (uo,Vo,ro) such that r which is the 
solution to (4), blows up in finite time, that is 

lim / r(x, t)dx — >■ oo, 

t->r*<oo Jq 

there exists a patch Hi C H, s.t for any single patch of measure greater than or 
equal to |f2! |, the modified model (8)-(9), with the same parameter set and initial 
condition, has globally existing solutions. In particular the solution r to (8), does 
not blow up in finite time. 

Proof. A proof of this conjecture can be established for certain special cases in one 
dimesion. Assume (4)-(6) blows up at time t = T* , for a certain parameter set 
and an initial condition r 0 (x ), such that r 0 ( n) = 0. We consider now introducing 
a refuge at the right end, starting at some positive a < tt. Now for x £ [a, tt] 
the density of v decreases, and the coefficient of r 2 is < 0, hence r is 

restricted on the right hand side of the domain. We assume that this is equivalent 
to introducing a Dirichlet boundary condition somewhere in the closed interval 
[a, tt]. Thus the modified model is equivalent to 

r 2 

(10) r t = d 3 Ar + cr 2 — w 3 -, on [0, 6], where a < b < tt, 

v + D 3 

and r x (0,t) = 0 and r(b,t) = 0. Here we assume the initial condition satisfies 
r 0 (6) = 0. Let us now compare (10) to 

(11) r t = d 3 Ar + cr 2 , on [0, b), where a <b < tt. 

The solution of the above is a supersolution to (10). However we know that there 
exists small data solutions to (11). This is easily seen via following the methods 
in [68]. Basically without loss of generality we may assume c = d 3 = 1. We then 
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multiply (11) by ry and integrate by parts in (0,6) to obtain 


d 

dt 




+ IMli = 0, 


thus 

We now define the functional 

^(*) = lkx(t)||I-|||r(t)||i, 

Since E'(t) < 0 we obtain 

^(i) < ,E7(0) = H^rolll - |||ro|||. 

Now we multiply (11) by r and integrate integrate by parts in (0, 6) to obtain 

(12) lj t \\ r\\l + \\ r x \\ 2 -\\ r\\l = 0 . 

This yields 

(13) ^ l|r|l ' + jB ^3 l|r|| 3 = 0 ' 

So If E{ 0) < 0, then E(t) < E( 0) < 0 and so 


(») 

This essentially yields 

(15) ^Iklll > I (lkllD f - 


Setting Y(t) := ||r(f)|| 2 - We derive the following differential inequality Y' > |y 3 / 2 . 
and so 


"Will > 


3||r 0 


3-t||ro||l' 


which means the solutions exist for t £ [0, X*), where 


(16) 


rj~i * 


3 

Ml’ 


r 0 ± 0. 


and then the solution r blows up at time T*. This means that if E( 0) > 0, then 
the initial data is actually small enough to ensure globally existing solutions [68]. 
What is required is 

(i7) IIMolli > §MII- 

This criteria can always be obtained for large enough refuge, that is, for 6 small 
enough. Since the boundary terms cancel, the norm here is equivalent to the 
I?q( 0,6) norm, hence by Sobolev embedding we have 

J \d x r 0 \ 2 dx >c(^J \r 0 \ 3 dx 



(18) 
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Now since f 0 \r 0 \ 3 dx << 1, for b chosen small enough we obtain 

(19) ' 

Thus combining (18) and (19) we obtain 


5 dx 


( 20 ) 


£ \d x r 0 \ 2 dx >c{j' |r 0 | 3 d^ > | (j* |r 0 | : 


dx 


This proves a particular case of conjecture 1. 


□ 


3.2. The Overcrowding Effect: Modified Model II. In high population den¬ 
sity areas a species should have greater dispersal in order to better assimilate 
available resources and avoid crowding effects such as increased intraspecific com¬ 
petition. These effects can be modeled via an overcrowding term, that has not 
been included in mathematical models of biological control. Consider the improved 
mathematical model that includes an overcrowding effect of the invasive species r, 
namely, 

r 2 

r t = d 3 r xx +cr 2 ~w 3 ——— + di{r 2 ) xx , 

v + D 3 

f uv \ f vr \ 

v t = d 2 v xx - a 2 v + wi ——— - w 2 [ — r . , 

\u + DiJ \v + D 2 ) 


( 21 ) 

( 22 ) 


with initial and boundary conditions as before. Again, the equation for u remains 
unchanged. The addition of di(r 2 ) xx = 4f( rr x ) x represents a severe penalty on 
local overcrowding. This is interpreted as movement from high towards low con¬ 
centrations of r, directly proportional to r. Hence, r attempts to avoid overcrowding 
and disperses toward lower concentrations. Such models have been under intense 
investigation recently and are referred to as cross-diffusion and self-diffusion sys¬ 
tems [53]. The mathematical analysis of such models is notoriously difficult [27, 28]. 
We limit ourselves to the one-dimensional case in the forthcoming analysis and its 
numerical approximations. Therefore the one dimensional Laplacian is considered 
in (6). 

The presence of blow up is not affected if we maintain Neumann boundary con¬ 
ditions at the boundaries. This can be seen in a straightforward manner. Let us 
assume the classical model, that is (4)-(6), blows up. Therefore, without loss of 
generality, there exists a <5 such that > <5 > 0. This implies that if we 

set H(t) = f n r(x,t)dx, then, 


d . . 6 

—Hit) > —= 
dt ~ ^\n\ 


H{tf 


leading to the blow up of Hit). Now, consider the integration over fl of (21). Since 
the overcrowding term integrates to zero then blow up still persists. However, a 
combination of a prey refuge with the overcrowding effect may prevent blow up. 
Further, we expect it takes a smaller refuge to accomplish the removal of blow up. 
This is precisely the following conjecture: 
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Conjecture 2. Consider the three species food chain model (4)-(6), a set of pa¬ 
rameters, such that c < j^, and an initial condition (no, t'o, 7'o) such that r which 
is the solution to (4) , blows up in finite time, that is 

lim / r(x, t)dx —> oo, T* < oo, 
t^T*<ooJ n 

there exists a patch ft 2 C £1, and an overcrowding coefficient cLp s.t for any single 
patch of measure greater than or equal to |S ^2 1 7 the modified model (21)-(22), with 
the same parameter set and initial condition, has globally existing solutions. In 
particular The solution r to (21), does not blow up in finite time. Furthermore, 
C , 2 C f2i, where fii is the patch found in conjecture (1). 

This is not difficult to see, as multiplying through by r and integrating by parts 
yields 



So If E( 0) < 0, then Eft) <15(0) <0 and so 

(23) ^IMli + J^r\r x \ 2 dx > ^||r||| > (||rj|i) 5 . 

Thus even if one has negative energy, E(t) < 15(0) < 0 so that (8) blows-up, due 
to the presence of the positive J n r\r x \ 2 dx term in (23), (21) will not blow up for 
small data, thus yielding a global solution. Thus there are initial data for which (8) 
blow ups but (21) does not, and so a smaller refuge would work in this case. The 
method follows by mimicing the steps in (12)-(20), with the additional J Q r\r x \ 2 dx 
term. 


3.3. Refuge, Overcrowding, and Role Reversal: Modified Model III. The 

prey refuge and overcrowding are included in the one-dimensional mathematical 
model. We also include a role-reversal of u within the protection zone of the refuge. 
This models the scenario in which two species may prey on each other in various 
regions where it is advantageous. Hence, the role-reversal of u will compete with the 
invasive species r. Figure 2 depicts the scenario of a one dimensional prey refuge 
for which outside the protection zone both u and r may predate on v. Hence, the 
model we propose of this scenario is given below, 


(24) r t = 
v t = 

(25) 

(26) u t = 


d 3 r xx T (c 


W 3 


(1 - b!(x))v + D 3 

d 2 V xx - a 2 v + h(x)wi ——— 
u + D i 


)r 2 + di{r 2 ) xx , 


-(1 - fol(x)) ( Wi 


vu 


v + D 2 


W 2 


vr 


v + D 4 


d\u xx + a\u — b 2 u 2 + (1 — bi(x))w 5 


vu 


v + A 


- b 1 (x)w 1 


uv 


u + D 3 


In light of this model we propose Conjecture 3: 


Conjecture 3. Consider the three species food chain model (4)-(6), a set of pa¬ 
rameters such that c < and an initial condition (uo,vo,ro) such that r which 
is the solution to (4), blows up in finite time, that is 

I r(x, t)dx —> 00 , T* < 00 , 
n 


lim 

t—>T* <00 
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FIGURE 2. This shows a plot of the refuge function bi(x) = 
|(l — tanh) 2 ^)). In certain regions u has a disadvantage over v, while 
outside of this u is able to prey on v. The invasive species effect on v is 
influenced by the prey refuge. 


there exists a patch H3 C fl such that the modified model (24)-(26), with the same 
parameter set and initial condition, has globally existing solutions. In particular, 
the solution r to (24) does not blow up infinite time. Furthermore, TI 3 C Tl 2 , where 
n 2 is the patch found in conjecture ( 2 ). 

All three conjectures are tested numerically in Section 5.3. 

4. Formal estimates and regularity 

4.1. Preliminaries. The nonnegativity of the solutions of (4)-(6) is preserved by 
application of standard results on invariant regions [57]. This is due to to the 
reaction terms being quasi-positive, that is, 

/ (0, v, r) > 0, g (u , 0, r) > 0, h (it, v, 0) > 0 , for all u,v,r > 0. 

Since the reaction terms are continuously differentiable on R +3 , then for any 
initial data in C (O) or L p (fi), p S (l,+oo), it is easy to directly check their 
Lipschitz continuity on bounded subsets of the domain of a fractional power of the 
operator J 3 (d\, d 2l dfif A, where I 3 the three dimensional identity matrix, A is the 
Laplacian operator and () denotes the transposition. Under these assumptions, 
the following local existence result is well known (see [17, 24, 45, 51, 57]). 

Proposition 1. The system (4)-(6) admits a unique, classical solution ( u, v,r ) on 
[0, T ln ax [x fl. If T m ax <C 00 then 

1 1™ 1 {IK*-OIL + IK*>OIL + IK*’OIL} = °°- 
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Lemma 4.1 (Uniform Gronwall Lemma). Let /3, C> and h be nonnegative functions 
in Lj oc [ 0, oo; R). Assume that /3 is absolutely continuous on (0, oo) and the following 
differential inequality is satisfied: 

^ < C P + h, for t > 0. 

If there exists a finite time t\ > 0 and some q > 0 such that 

rt-\-q /*£+<? 

/ C (r)dT < A, / p{r)dr < B, and / h(r)dT < C, 

Jt Jt Jt 

for any t > t\, where A , B, and C are some positive constants, then 


m< 


B 


+ C I e A , for any t > t 3 + q. 


4.2. Improvement of Regularity of the Classical Model. 


4.2.1. Improvement of the Global Existence Condition. In this section we improve 
the global existence conditions that has been derived recently in [43]. Here, we 
provide the previous result, Theorem 4.2. 

Theorem 4.2. Consider the three-species food-chain model described by (4)-(6). 
For any initial data (uo,vo,ro) € L 2 (fl), such that 11 vo 11 oo < — D 3 , and param¬ 
eters (c, W 3 ,D 3 ) such that ||r||oo < — D 3 , there exists a global classical solution 

(it, v, r) to the system. 

We now give an improved result for the ODE case, which is easily modified 
for PDE case. In essence, given any initial data (however large), there is global 
solution, for 02 appropriately large. This is summarised in following lemma, 

Lemma 4.3. Consider the three-species food-chain model described by (l)-(3). As 
long as c < given any initial data (uo,vo,ro), however large, there exists a 
global solution ( u,v,r ) to the system as long as the parameter <22 is s.t 

a,> Wl+ clr 0 \ln( -]fL_ y 


Proof. Consider the following subsystem 


(27) 

dv 

— = —02V + w\V 

dt 

(28) 

dr 9 

— = cr . 

dt 


The v here grows faster than the v determined by (2). Similarly, the r here blows 
up faster than r in (1). In order to drive v down below w^/c — D 3 , we use the exact 
solution to (27), that is, 

v = e- {a 2 ~ Wl)t v 0 

Assuming that 02 > W \, then for v to be below ir 3 /c — D 3 implies that 

e -( 02 - Wl ) t Vo < W3 / c _ D3 

Thus when t > T* = — 1 — In ( —) j w e shall have 

0.2 —Wi \W3/C—D 3 J' 

e -(o 2 -Wl)t 


Vq = V < w 3 /c - D 3 . 
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Now (28) blows up at time T** = . So if we choose 02 appropriately, we can 

make T* < T** = d. This is done by choosing 


1 


In 


M 


a 2 — wi \w 3 /c - D 3 


= T* <T* 


1 

cM 


Thus if 02 is s.t 02 > w\ + c|ro| In Ds J , then v < w 3 /c — D 3l before r blows 

up. Therefore, if we consider the r determined by (1) then this will certainly not 
blown up by time t = T*. At this point however, c — is negative, so r in (1), 

cannot blow up after T*, and will decay. In this latter case the global attractor is 
a very simple, (u,v,r) = (u*,v*, 0). Hence, we have an extinction state for r. 

In short, we see that for any initial condition that is specified, there exists an 02 
(depending on the initial condition), s.t (l)-(3) has a global solution. □ 


4.2.2. Uniform L 2 {Vi) and H 1 (12) Estimates. We recap the following estimate from 
[42] 

CC X 


MI2 < e c i t ||u 


Therefore, there exists time t\ given explicitly by 


61 

such that, for all t>ti, the following estimate holds uniformly: 


Mia < 1 + 


CC 1 


Here t± is the compactification time of ||m|| 2 - 

We also recall the following local in time integral estimate for V«, from [ 12] 

rt +1 1 1 rt +1 1 / CC 1 \ C 1 

l nvuf 2ds <- M mi +T J t Cts <-{ i+ ^) + -, lor t> h . 

Now we move to the L 2 estimate for the v component. This cannot be obtained 
directly. Thus we use the grouping method again by multiplying ( 6 ) by w 1 and (5) 
by wq, adding the two together, and setting w = W\u + wqv , we obtain: 

(29) Wt = + (di — d 2 )wiAu + w\a\u — w\biu 2 — woa 2 V — W 0 W 2 ' 

We add a convenient zero, —a 2 W\u + a 2 W\u, to (29) to obtain 

wt = d 2 Aw + ( d± — d 2 )w\Au + wia 3 u — w\b\u 2 — a 2 W + a 2 W\u — W 0 W 2 


v + D 2 


v + D 2 


By multiplying (29) by w, integrating by parts over 12, and keeping in mind the 
positivity of the solutions, we find that 

\j f \\M\l+ d 2\Ww\\l+ a 2\\w\\l 

< wi(ai + 02 ) / uwdx + / {d 2 — di)wiVu ■ Vwdx. 

Jft JQ. 
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This yields 


— ||w||§ + 2min(a 2 ,d 2 ) (||Vw;|| + |M| 2 ) 


< 2wi(cii + 02 ) / uwdx + 2(d 2 — d x )w x / X/u-Vwdx 

J Q J Q 


< min(a 2 , <^ 2 ) I |tc|I 2 + 


2 (2wi(ai + a 2 )) 2 ||2 . , ... ||2 


2 min(a 2 , d 2 ) 


|u || 2 + min(a 2 , d 2 )||Vu;| 


such that 


+ ( 2 ( <i 2 -d i y 

2 mm(a 2 , a 2 ) 

d 


dt 


w Hi + min(o 2 , d 2 ) (||w || 2 + 11 Vic|||) 


2 min(a 2 ,d 2 ) 
Now via estimates in [42] we obtain 


< ( 2 "'i(-+-o) 3 ii„|g+ (y a -. ii ‘)n,) 2 ii?.ii. 


2 min(a 2 ,d 2 ) 


Hl2<C, 


|o, In f — 


Molll+molkolll 


C 


)} 


for t > £3 where £3 is explicitly derived in [42] , t 3 = max 
which implies that 

IMI 2 < C, for t > t 3 . 

Notice by multiplying (5) by v and integrating by parts one obtains, 

\j t \\ v \\l + d ' 2 \\^ v \\l + a 2 \\v\\l < w x \\v\\l 
Integrating the above on [t 3 ,t 3 + 1] yields 

/ ||Vu|||ds< / w 1 \\v\\lds+ -||w(f 3 )|| 2 < (w x + 1)C < C 

Jt 3 Jt 3 1 

Thus via the mean value theorem for integrals we obtain a time t 3 € [t 3 ,t 3 + 1] s.t 

\\Vv(t* 3 )\\ 2 2 ds<C 

In addition, we recap the uniform H 1 estimates found in [42], 

(30) \\Vu\\l<C, 


for t > Ti, the H 1 compactification time of the solution. We next need to derive 
H 1 estimates for the v component. This is tricky, since when we multiply (5) by 
—Aw, integrate by parts, and then apply Young’s inequality with epsilon we obtain 

^l|Vw || 2 + a 2 \\Vv\\l < wi\\v\\l + w 2 \\r\\l 

However, there is no global estimate on r, as it may possibly blow up, but we know 
that there is always a local solution on [0, \T*] : where T* is the blow up time 
of r t = d 3 Ar + cr 2 . Furthermore, on this time interval the solution is classical. 
Thus we make local in time estimates on this interval, using standard local in time 







16 


BEAUREGARD, BLACK, PARSHAD, QUANSAH 


estimates on IHloo from [48]. Using Gronwalls inequality via integration on the 
time interval [t 3 ,t] yields, 

|[Vu||l < e- 2 a ^\Vv(tml+-(wiC + w 2 C 1 ^^) 

< 1 + — (w 1 C + w 2 C 1 ^^ 

a 2 \ d 3 

for t>T A = max jo, ; |. 

Note, Lemma 4.3 requires us to manipulate a 2 , in order to prove large data global 
existence. From a practical point of view this is not always easy as manipulating 
the death rate directly may not be condusive to realistic control strategies. We 
next provide the following improved result that is valid even if ||uo||oo > -yf — D 3 , 
without manipulating a 2 . 

Theorem 4.4. Consider the three-species food-chain model described by (4)-(6). 
For certain initial data (uo,vo,ro) £ L 2 (Cl), there exists a global classical solution 
(u,v,r) to the system, even if Huolloo > qy — D 3 , as long as the parameters are 
such that ||u| |oo < -yf — D 3 , by the H 2 absorption time of v ( denoted Tq), If Tq 
satisfies Tq <T*, where T* is the blow up time of the following PDE 

dtr — d 3 Ar = cr 2 . 

with the same initial and boundary conditions as (4) 

In order to show that blow up can be avoided as stated above, even in R 2 our H 1 
estimates via (30) must be improved. This is because the embedding of H 1 L°° 
is lost in R 2 . Thus, if we are to control the L°° norm in R 2 we need H 2 control. 
This is achieved in the next subsection. 

4.2.3. Uniform H 2 (Ui) Estimates. We will estimate the H 2 norms via the following 
procedure, we rewrite (6) as 

J A , 9 ( uv 

u t — d\Au = aizt — b 2 u" — w 3 - 

\u + D o 

We square both sides of the equation and integrate by parts over f l to obtain 

ll u t||i + di11 Aw11 2 + 11 Vu| 1 2 

< C ((ai) 2 ||u || 2 + (zc 0 ) 2 11p11 2 ) + (^ 2 ) 2 11^|I 4 < C. 

This result follows by the embedding of H 1 4 L 4 L 2 . Therefore, we obtain 

IKII 2 + ^l||^ u ll2 + ^l|Vu||2 < U. 

Due to classical local in time regularity results, see 1, the solutions are C' 1 (0,T; C 2 (U))), 
thus 'j f 11Vzz.1 12 cannot escape to — 00 , and is bounded below, upto the existence time. 
Thus we obtain 

11Azz|11 < C 2 (t). 

These constants may depend on t since even though ^||Vu|| 2 cannot escape to —00 
in finite time the case of it decaying like —t is not precluded. We now make uniform 
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in time estimates of the higher order terms. Integrating the estimates of (32) in 
the time interval [T,T + 1], for T > T 5 = max(T 1 , T 4 ), where T 5 is the maximum 
of the H 1 compactification times of u, v yields 


pT-\-l 

J \\u t \\ldt < Ci, 



11 Ait] < Ci- 


Similarly we adopt the same procedure for v to obtain 


pT~\~ 1 

J \\v t \\ldt < Ci, 



||A«||ldt<C 2 . 


However, the above estimates for v are made on the time interval [0, A-]. 

Next, consider the gradient of ( 6 ). Following the same technique as in deriving 
(32) we obtain for the left hand side 


11Viif1 12 + di11V(Au)|I 2 + -j- 11Ait]1 2 + f AuVufndS 

dt " Jon 

= W^utWl + d\ 11V(An)11 2 + ^\\Au\\j + A u|(Vu • n )dS 

(31) =\\Vu t \\* + d 1 \\V(Au)\\t + j t \\Au \\ 2 2 
This follows via the boundary condition. Thus we have 

MVittlll +d 1 ||V(A«)||5 + 

= ( v (ain-M 2 -R>o(^))) 

(32) < C'(||u||| + IMI 4 + |]Au|I 2 + 11Av|I 2 ) 


This follows via Young’s inequality with epsilon, as well as the embedding of 
H 2 (fl) c —>• bF 1 , 4 (fl). This implies that 

— ||Au||2 < C|IAit|I 2 + C'dlulll + |M|t + 11 At7| ||) 

We now use the uniform Gronwall lemma with 


m = ha«iii, m = c, h{t) = c(\\u\\t + nun 4 + im 2 ), q = 1 

to obtain 

11Art11 2 < C, ||Au || 2 < C, for t>T 6 = T 5 + 1 

The estimate for v is derived similarly, but again we assume we are on [0, ^-]. Thus, 
via elliptic regularity, 

(33) IMIit 2 < C, |M|lp < C, for t > T e 

Since H 2 c —> L°° in R 2 and R 3 , the following estimate is valid in R 2 and R 3 

(34) Halloo < C, ||u||oo < C, for t > T 6 

Thus even if the data is such that 11 w 0 11 00 > yy — D 3 , if the parameters and the 
data (r 0 ,u 0 ,Uo) are such that — — D 3 < C 3 C, (for the C in (34), where C\ is the 
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embedding constant for H 2 [fd) c —>■ L°°{ f2)) and 


T 6 (r 0 ,v 0 ,u 0 ) < l + max<jO,Ti, 

d 3 


ln(||Vi;(^)|||) /willuolll + w 0 IL '- 112 


, In 


a 2 


C 


+ 1 


< 


< T , 


then we obtain 

IMloo < — - £>3 < c, for t > T 6 . 
c 

Since T e < T* the solution to d t r — d 3 Ar = cr 2 , has not blown up yet, and 
since this is a supersolution to r solving (4), this r has definitely not blown up. 
However, from this point on the coefficient of r 2 in (4) is negative, as — v +b 3 ) < - 

— < 0, hence r solving (4) can never blow up, from this point on. 

This proves Theorem 4.4. 


4.3. Long Time Dynamics. 


4.3.1. Gradient Estimates for the Time Derivatives. Now we assume there are glob¬ 
ally existing solutions, and we aim to investigate the regularity of the omega limit 
set. Hence we assume we are working with initial data and parameters for which 
r is globally bounded. To begin, we obtain useful integral in time estimates. In 
particular, we integrate (32) in time from [T, T+l] where T > Tq, the H 2 absorbing 
time, to obtain 
/•T+l 

j T \\Vu t \\ldt < C(\\u\\i + IMI2 + ||Ait || 2 + ||Au|| 2 ) + ||Au(T )|| 2 < C. 

Similarly, an estimate for Vu* can be established, namely 



||V« t || \dt<C. 


We now develop higher order estimates for the time derivatives. First, consider 
the partial derivative w.r.t t of equations ( 6 ), multipling by —Au t , and integrate 
over fl we obtain 


d 

dt 


11V It* 11 § 


+ d\ 11 Aut 115 


o-i 11V rt* 111 + b 2 


/ u(ut)Autdx + wq 
J n 


+w 0 / v(u t )Au t 


Dn 


(u + Dq) 


:dx. 


- —v t Au t dx 

u + D o 


Using Holder’s inequality, Young’s inequality with epsilon, and our earlier estimates 
yield 


^l|Vu t || 2 + di||Au t || 2 


< 


ai\\Vu t \\- 2 + C'Hull^llutlla + -^||Au t || 2 + C'i I I ll 

+-^\\A u t\\l + Cal |w| |oo I I'Ut I li + -^\\Au t \\l. 


Thus we obtain 

^;||Vu t ||2 < ai||Vu t ||2 + C'||M||^ 0 ||u i ||2 + UillUfll; + C'2||'y||oo||ut|l2- 
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The application of the uniform Gronwall Lemma with 
0 (t) = ||VittHl, C(t) = ai, 

h(t) = C\\u\\ 00 \\u t \\l + Cx\\v t \\l+C 2 \\v\\ 00 \\u t \\l, q = 1 
gives us the following uniform bound 
(35) ||Vut||2< C, t >T 7 = T 6 + 1 

Similar methods applied to the equation for v and r yield 

llVwtH 2 2 <C ,t>T r , ||Vr t | \\<C, t>T r . 


4.3.2. Regularity of the omega limit set. We now show that the asymptotic state 
(0, v , u) under certain parameter restrictions is an attractor with more regularity 
than was derived previously. Essentially, we consider (ro,uo,uo) £ L 2 {Pt) for which 
we have globally existing solutions, that is the maximal existence time for the 
solutions T = oo. Then we consider the omega limit set for such solutions. 

(36) w(r 0 , v 0 , u 0 ) = (J {(r(s), v(s), u(s)) : s > t} 

t> o 

We state the following result. 

Theorem 4.5. Consider the three-species food-chain model described by (4)-(6). 
For initial conditions (ro, Vo, u o), such that there is a globally existing solution, there 
exists a parameter a 2 , (depending on the initial condition), for which the omega limit 
set uj(ro,vo,uo) = (0 ,v,u). Furthermore this is an attractor with H 2 (Q) regularity. 

Remark 5. Thus for suitably chosen < 22 , the attractor for such solutions is the 
extinction state of r. 


Proof. We have shown that the system is well posed, under certain parameter, 
and data restrictions, via theorem 4.4. Thus, there exists a well-defined semigroup 
{S(t)} t>0 : H —► H. Also, via Lemma 4.3 we can find an 02 such that r decays to 
zero. The estimates via (33) establish the existence of bounded absorbing sets in 
H 2 . Now we rewrite the equation for a sequence u n as follows 

diAu n = d t u n - a±u n + b 2 u 2 n + w 0 ( ) , 

\u n + D 0 J 

We aim to show the convergence of both the right hand side strongly in L 2 . Due to 
the uniform bounds via (30), (35), (33) and the embedding of R 1 (D) c —> i 4 (ff) c —^ 
L 2 (Q) and H 2 (Q) c —> L°°( Cl) we obtain a subsequence u nj stil labeled u n s.t 

d t u n —A d t u*, 

L 2 

ai u n —> aiii*, 
b 2 ul ^ b 2 (u*) 2 


as L \ (u n — u*)(u n + u*)\ 2 dx < llun+u* 


\u n -u 


0 , as u n 


u*. Similarly, 


Wo 


f U n V n \ 
\Un + Dq J 


w 0 


f u*v* \ 

Vu* + aJ 
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L 2 

as v n —> v*. Thus, given a sequence {'Wn(0)}^L 1 that is bounded in L 2 (fl), we 
know that, for t > T 7 , 

S(t n )(Au n (0)) —► Au* in L 2 (fl). 

This yields the precompactness in H 2 { 11), hence the closure in (36) is in H 2 (fl). 
A similar analysis can be done for the v and r components. Hence the theorem is 
established. □ 


5. Numerical Approximation 

In this section a numerical approximation for the one and two dimensional mod¬ 
els are developed in order to numerically demonstrate our earlier conjectures ( 1 ), 
(2), and (3) in addition to exploring the rich dynamics the mathematical models ex¬ 
hibit. A one-dimensional spectral-collocation method is developed to approximate 
the one dimensional equations (24)-(26) which include overcrowding, refuges, and 
role-reversal. We then offer a second order finite difference approximation of the 
two-dimensional equations ( 8 ), (9), and ( 6 ) that investigate the effect of a refuge. 
Overcrowding and role-reversal effects are not investigated in the latter numerical 
approximation. The two dimensional approximation employs a Peaceman-Rachford 
operator splitting and techniques developed by one of the authors in [5]. This ap¬ 
proach takes advantage of the sparsity and structure of the underlying matrices and 
the known computational efficiency of operator splitting methods. Without loss of 
generality the domain is scaled and translated to (— 1 , 1 ) and (— 1 , 1 ) x (— 1 , 1 ) in 
one and two-dimensions, respectively. 

5.1. Approximation in One Dimension. We develop a spectral-collocation ap¬ 
proximation of (24)-(26). These equations can be written compactly as 


w t = £w, 


where w = ( u,v,r) T and C is the differential operator that includes the Laplacian 
operator and reactive terms. 

The spatial approximation is constructed from a Chebychev collocation scheme 
[13, 25, 52]. The spatial approximation is constructed as a linear combination of 
the interpolating splines on the Gauss-Lobatto quadrature. The resulting system 
is then integrated in time using an implicit scheme, in particular a second order 
Adams-Moulton method. The Chebychev collocation approximation allows for a 
high order spatial approximation [25]. This offers the ability to capture fine res¬ 
olution details with a relatively smaller number of degrees of freedom. However, 
there is a downside. The resulting matrices are full and ill conditioned [52] . This is 
problematic for the inversion in the resulting linear solve. Nevertheless, the spec¬ 
trum associated with the second derivative operator is real, negative, and grows in 
magnitude like 0(N 4 ) [23, 25]. 

The collocation scheme is constructed on the Gauss-Lobatto abscissa with re¬ 
spect to the Chebychev weight, 

(37) 
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where j = 0,..., N. An approximation is then constructed as a linear combination 
of the Lagrange interpolates on the abscissa, 

N 

(38) wjv {x,t) = a 

j=o 

where wn(x, t ) is an approximation to the unknown populations and SLj(t) is a 3 x 1 
vector of fourier coefficients. The basis functions are defined by 




(_l)Ar+j+i(l -x 2 )T' n (x) 
CjN 2 (x - xj ) 
f 2 if j = 0, N, 

] 1 otherwise. ’ 


and each <f>j(x) is a polynomial of degree N [25]. The basis functions also satisfy 
4>j{xi) = Sji. This is used to develop an approximation to our differential operator 
and a system of equations is constructed via a discrete inner product, 


<W t ,0 m >AT = < P N Cw N , (jim > N , 


where Pn is the projection operator onto the space of polynomials of degree N. The 
discrete inner product is given below, and it is based on an inner product defined 

by 


(39) 


<u,v> w 



u(x)v(x) 

\/l — X 2 


dx. 


The integral can be approximated as a sum over the Gauss-Lobatto quadrature for 
the Chebychev weight and is exact for polynomials up to degree 2 N — 1, 



p 2 N-i(x) 

a/ 1 — X 2 


dx 


N 

y^ J P2N-l{x k )w k . 

i =0 


The resulting norm results in an equivalent norm for polynomials up to degree 2N 
[8]. The abscissa, Xi, are the same as those defined in equation (37), and the Wi are 
the quadrature weights. The resulting numerical approximation is then constructed 
using this finite inner product 


N 

< f,g>N = ^2 f( x i)g( x i) w i- 

i =0 


The resulting approximation is an extension of that given in (39) and is constructed 
using a Galerkin approach, 


< w t , <j> m >n 


< PNjCwNj^m >N, 


E 


d_ 

dt 


WA r(Xi,t)<l>m(Xi)Wi 


N 

E Pjv£ WAr,^m(a; i)Wi. 


2 — 0 2=0 

The basis functions, r), are Lagrange interpolants, and after substitution of the 
definition in equation (38) the result can be greatly simplified, 


a m (0 = p n£wn 


where m = 0,..., N. 


X—X. 
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The Adams-Moulton method used to advance the initial value problem in time 
necessitates solving a nonlinear system of equations at each time step. This is 
achieved through a Newton-Raphson method. The ensuing linearized problem in 
the Newton-Raphson method is solved through a restarted GMRES [37]. 

5.2. Approximation in Two Dimensions. A second order approximation of 
(8), (9), and (6) is developed. While the spectral-Galerkin appoximation of the 
previous section may be extended to two dimensions, the computational cost of the 
nonlinear solve is troublesome due to the sensitivity of the matrices and the lack of 
their sparsity. Here, the approximation is still of high order, while the underlying 
matrices are block tridiagonal or tridiagonal with diagonal blocks. This enables 
direct inversion techniques, in particular the Thomas algorithm. 

The governing equations are written compactly as 

Wf = D Aw + f, 

where w = (it, u,r) T , f= (/(w,x, t), g(w, x,t), h(w,x,t)) T are the reactive terms 
that depend on space, time, and the species it, v, and r, x = (x,y), Aw is taken 
component-wise, and D is a diagonal matrix with nonnegative entries g?i, and 
c? 3 . From semigroup theory the formal solution is 

w(x, t) = exp (tDA) w(x, 0) + [ exp {(t — r)A} f(x, r)dr, 

Jo 

where exp{} is the evolution operator associated with D. A suitable quadrature 
is used to approximate the integral. Here, a second order trapezoidal rule is used, 
that is, 

w(x, t) = exp (t.DA) ^w(x, 0) + |f(x, 0)^ + |f(x, t ) + 0(St 2 ). 

This motivates the implicit method, 

w(x, tk+i) = exp(tDA) ^w(x,t fc ) + ^f(x,t fc ) 

+ yf(x,4+i) + 0(6t 2 ), 

where tk+i = tk + St. The exponential is approximated through a Peaceman- 
Rachford operator. This creates the second order implicit method, 

w (x, tk+i) = {I- ^A)~ 1 {I+ ^A){I+ ^j-B) (w(x,4) 

+ ^ f ( x ^fc)^) + y f ( x ,^+i) + 0(St 2 k ), 

o2 o2 

where Stk is the variable time step, A = -§^2, and B = . The last reactive term 

does require the solution at the next time step. To avoid a required nonlinear solve, 
we use a first order approximation to avoid this. This simplification maintains the 
second order accuracy of the approximation. We write the solution in an alternative 
form, 

(7-^A)(J-^H)w(x,4 + i) = (I+ S yA)(I+ S yB)(w(x,t k )+ 6 ^{(x,t k ^j 

■+—~ 1 ) + 0(St 2 k ), 
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This can be conveniently solved through ADI procedures [59], that is, by splitting 
the problem into, 

{I - ^A)w(x, t k+ 1) = (/ + ^B) w(x, t k ) + ^f(x, t k ) 

(/- ^B)w(x,t fe+ i) = (/ + ^A)w(x,4) + < y^f(x,4 + i). 

We see that the first equation keeps A implicit while B is explicit. We then take 
our intermediate solution, w, and solve the second equation keeping B implicit 
and B explicit. Now, the spatial operators may be approximated through the 
two-dimensional Chebyshev approximation similar to that of the previous section, 
however we choose second order central differences. Let Xi = — 1 + ih x and i/j = 
— 1 + jh v , where h x = 2/(N — 1), h y = 2/(M — 1), for i = 0,..., N — 1, and 

j = 0, ...,M — 1. Let Ah and B k be the second order approximations to the 

operators A and B. The approximation is utilized throughout the entire two- 
dimensional domain. At the boundary, we require an equation for w at ghost 
points, that is locations of X-i, y~i, Xm +i and i/M+i- These are established using 
central difference approximations to the Neumann boundary conditions. We then 
substitute these back into the system of equations. This maintains the second 
order accuracy and the tridiagonal structure of the equations. At each step, the 
tridiagonal equations in the ADI method are directly solved using the Thomas 
algorithm which comes at a expense of 0 (NM). 


5.3. Numerical Results. We numerically demonstrate that there is good evidence 
for conjectures (1), (2), and (3). In the one dimensional setting, we translate and 
scale the domain to (0,7r) and use a refuge function of 


h(x) 


1 - tanh (%f) 

2 


The refuge size is delineated by the value a since that is the location where the 
gradient is steepest. As expected the blow up time is affected by the presence of a 
refuge. For instance, for a fixed parameter set we look at the effect on the blow up 
time with a refuge, refuge with role reversal, and role reversal and overcrowding as 
compared to the classical model. The results are shown in Figure 3. 

Interestingly, in the modified model with role reversal, it is found that blow up 
does not occur when b\(x) = 0 or b\{x) = l,Vx in the spatial domain. 

We also see that the blow up times are influenced by the size of the refuge. In 
fact, there is critical refuge size for which given any refuge size greater then the 
solutions will not blow up. This is evidenced in Figure 3 by the steep gradient 
of the curves around 2.8. Therefore, the population of the invasive species can be 
controlled. In situations where there is only a spatial refuge the blow up time curve 
is always increasing. This is consistent with our two dimensional results. 

It is clear that there exists a critical refuge size for which refuges larger than 
this will prevent blow up. In the one-dimensional model, without role reversal 
or overcrowding effects, we investigate the effect of the critical refuge size versus 
the size of the initial condition of v while maintaining the other initial conditions 
and parameters values. We use the same parameters as given in Figure 3 and 
vary the uniform initial condition on v. Interestingly, Figure 4 shows a logarithmic 
dependency of the critical refuge size to prevent blow-up of r, on the initial condition 
size of v. 
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FIGURE 3. This shows blow up times for the classical model ( 61 ( 2 :) = 0 
and d ,4 = 0 ) versus the modified model with the various biological con¬ 
trols. One sees that there is a critical patch size, such that for patches of 
length greater than this, there is no blow up. Inclusion of the overcrowd¬ 
ing term, greatly decreases this critical size. The parameters used are 
ui =■ 1, 0,2 — 1, — 0.5; Do = 10, D\ — 13, D 2 — 10, Do — 20, D 4 = 

D 2 , c = 0.055, wo = 0.55, wi = 0.1, W 2 = 0.25, W 3 = 1.2, W 4 = 
100, W 5 = 0.55, do = do = d\ = .1, do = k\c — |j||, where k is given 
in the plot. The initial conditions used are u(x, 0) = r(x, 0) = 10 and 
u(a:, 0 ) = 2000 . 128 grid points were used with a temporal step size of 
10 -3 . In the case of overcrowding and role reversal, as we increase the 
value of k similar results were observed. 


The size of the spatial refuge has an influence on the blow up time. However, 
the location and number of spatial refuges also influences the blow up time. In fact, 
blow up is sometimes eliminated depending on the initial condition, refuges, and 
their subsequent locations. For instance, if we consider a uniform initial condition 
in r and compare the blow up time in a situation of a single refuge of width located 
near the boundary versus the case of evenly splitting this refuge then we find the 
blow up time is increased. This delay is exasperated if the two refuges are farther 
apart. Of course, if we consider a different initial condition then blow may not 
just be delayed, rather removed! For instance if we consider a concentrated initial 
population of r for which the highest concentration of r contains the spatial refuge 
then blow up can be eliminated! Hence, the concentration of r inside the refuge 
may protect the species enough so that the population of r decays sufficiently to 
avoid blow up in its population. 
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FIGURE 4. This shows a plot of the critical patch size for a given 
initial uniform condition size of v for which blow up in the invasive 
specie population r will occur. Identical parameters and resolution are 
used as in Figure 3. The domain was translated and scaled to ( 0 , 7 r). 


To illustrate this consider the two-dimensional model ( 8 ), (9), and ( 6 ) with an 
initial condition of 


u(x,y, 0) = cos(27ra) cos(27r?/) + 30, 

v(x, y, 0 ) = u{x, y, 0 ) + 200 , 

r(x,y, 0 ) = 100 exp (— 10 (x 2 + y 2 )), 


with parameters d\ = c ?2 = do = .1, ai = 5, 02 = .75, 62 = -5, wq = .55, w\ = 1, 
W2 = -25, UI3 = 1.2, c = .055, Do = 20, D \ = 13, D2 = 10, D3 = 20. Clearly, the 
coefficient c— < 0. If bi(x,y) is zero throughout the entire spatial domain there 
exists blow up in the r population. However, if we have a circular refuge such that, 


h{x,y) 


1 x 2 + y 2 < R 
0 else 


for R = .5 then blow up is avoided. Figure 5 shows the total population versus time. 
We can see the population starts to increase rapidly, but the increase is attenuated 
as a result of the spatial refuge. Hence, the population growth is not sustained 
and begins to decrease. Hence, the size and location of the refuge and the initial 
condition play a delicate balance in preventing blow up. 

In two dimensions if we choose a particular shape of the spatial refuge it was 
conjectured that there exists a critical size for which blow up is avoided. Here, we 
look at two situations: a square and circular refuges with increasing size centered 
in the middle of the spatial domain. We let fri(x) = 1 inside the refuge while 0 
outside. This clearly delineates the protection zones. For a fixed parameter regime 
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Time 


Figure 5. The population of r(x,y,t), that is, f n r(x,y,t)d$l, is 
shown versus time. With a circular refuge, we see that blow up 
is avoided. This provides experimental evidence that the location 
and size of the refuge is important to avoid blow up in the invasive 
population. The dashed line represents the population in the case 
where there is no spatial refuge. The simulations were done on a 
50 x 50 uniform grid with a temporal step of .0001. 


and initial conditions of 

u(x, y , 0) = r(x, y , 0) = cos(27rx) cos(27ry) + 30, 
v(x, y, 0) = cos(27rx) cos(27 ry) + 230, 

we determine the critical refuge size. Each calculation is carried out with a 50 x 50 
grid. The temporal step was fixed at .001. The parameters used: d\ = cfe = do = -1, 
a± = 1, b\ = .5, wq = -55, w\ = 2, = .25, W 3 = 1.2, c = .055, Dq = 20, Di = 13, 

D 2 = 10, D3 = 20. It is found that for a square the critical refuge area is 2.5004, 
roughly 62.5% of the spatial domain. The critical refuge area for the circular refuge 
is approximately 2.6661, roughly 66.7% of the spatial domain. 

6. Spatio-temporal Dynamics 

6.1. Turing Instability: Effect of Overcrowding on Turing Patterns. In 

this section we shall investigate the effects of overcrowding in the absence of refuges 
and role reversal, in the classical model. Therefore, we focus on whether an appro¬ 
priate choice of d 4 can induce Turing Instabilities. It is shown in [43] that diffusion 
processes can destabilize the homogenous steady state solution when d 4 = 0. For 
convenience, we restate the one-dimensional model given in equations (6), (21), and 
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( 22 ), 

du 

~dt 


(40) 

= d\U xx + a\u — b 2 u' 

(41) 

dv 

dt 

= d 2 v xx - a 2 v + W! ^ 

(42) 

dr 

dt 

= d 3 r xx + cr 2 w 3 

V 


UV 


UV 


U + D i 
‘2 


U + Dq 
~ W 2 


V + D 2 


d^{v ) X X‘ 

Consider the linearization of (40)-(42) about the positive interior equilibrium point 


= ( u*,v*,r* 


where u* = ai 2 ^'° 0 + \j ( ai 2b 2 -, D °) 2 ~ (~ ° t '* 6 , aiDn ')> v * ' s the the spatially homoge¬ 
nous steady state solution, and r* = v ( ^ 1 + Di — a 2 ) (see [43]). For instance, 
for the parameters given at the beginning of this section u* = 10.110031, v* = 10, 
and r* = 2.997897. Consider a small space time perturbation, that is, 


W = U-U* = O(e), where e -> 0, 

with U* as the positive interior equilibrium point given as E§ and U = ( u,v,r ). 
Substituting and collecting linear terms of order O(W), we obtain 

dW 


where 


dt 

AW, • n = 


D, 


DAW + JW, 

0 for x € dfi, i = 1, 2, 3. 


di 

0 

0 


0 

d 2 

0 


0 

0 

o ?3 + 2d^r* 


is the diffusion matrix and 
— b 9 —{— 


(u*+D 0 y 

v* Diwi 

( u *+ D !)2 

0 


u* -\-D q 

v * W 2 r* 

(■ V*+D 2 )2 
‘ 2 w a 


v*+D 2 

0 


An 

A 12 

CO 

a 2 i 

^22 

a 23 

CO 

^32 

a 33 


(u*+-D 3 ) 2 

is the Jacobian matrix associated with the ordinary differential equation part of 
model (40)-(42). 

Let 

wo 


W (e,i) = 


\t-\-ike 


where e is the spatial coordinate in f2, Ui(i = 0,1,2) is the amplitude, A is the 
eigenvalues associated with the interior equilibrium point, Eq and k is the wave 
number of the solution. Upon substituting, we obtain the characteristic equation 

(43) |J- AI-/c 2 D| =0, 


where I is a 3 x 3 identity matrix. The sign of Re( A) indicates the stability, or lack 
thereof, of the equilibrium point Eq. The dispersion relation is 

P( A) = A 3 {k 2 ) A 3 + A 2 (k 2 ) A 2 + A ± {k 2 ) A + A 0 (k 2 ). 
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The coefficients of P( A) are determined by expanding (43), namely, 

A 3 (k 2 ) =1, 

A 2 (k 2 ) =(di + d 2 + {d 3 + 2d^r*) )k~ — Hu — H22 — H33, 

Hi(fc 2 ) —^22^33 — ^22(^3 + ‘2d^r*)k 2 — H23H32 — 4 li 2 ^ 42 i — d 2 A 33 k 2 

+ d 2 {d 3 + 2 d 4 r*)fc 4 + f (difc 2 - H n ) (d 2 k 2 + (d 3 + 2 d 4 r*)k 2 



Ag(k 2 ) — ^ [dik 2 — An) (H22H33 — 41.22(^3 + 2 d 4 r*)k 2 — H23H32 — A 33 d 2 k 2 
+ d 2 (d 3 + 2 d 4 r*)k 4 )'\ — Hi2H2i(d3 + 2 d 4 r*)k 2 + Ai2j421^33- 


To check for stability of the equilibrium solution Eg, we use the Routh Hurwitz 
criterion. This state that for Eg to be stable we need 

(44) A n (k 2 ) > 0, Vn and Ai(k 2 )A 2 (k 2 ) > H 0 (fc 2 ). 

Contradicting either of these statements ensures instability for Eg. Finally, for 
diffusion to cause a Turing instability it is sufficient to require that around the 
equilibrium point we have 

Re(A (k = 0 )) < 0 , and 

Re(A(fc > 0)) > 0. 

We refer the reader to a well detailed analysis of this in [22]. Hence, for a Turing 
instability to occur, we require that (44) is satisfied when k = 0 (without diffusion) 
and at least one of the equations in (44) changes sign when k > 0 (with diffusion). 
By this we consider letting Ag(k 2 ) become negative when k > 0 and satisfy (44) 
when k = 0. This suggests that spatial patterns should be observed. Our parameter 
search has not yielded a parameter set for which (44) holds while d 4 = 0. When 
d 4 > 0, at least one of these inequalities changes sign. Thus we cannot conclusively 
say that d 4 will induce or inhibit Turing pattern. However, it is conclusive that 
d 4 certainly has an effect on the type of patterns that do form. In particular, the 
patterns fall into two types: spatial patterns and spatio-temporal patterns. The 
conditions for which type forms is succinctly described via Table 1. 


Case 

Aq 

A 2 Ai — Ho 

Hi 

Pattern Type 

1 

+ 



spatio-temporal 

2 


+ 

+ 

fixed spatial 


Table 1. Conditions on the signs of coefficients for different types 
of patterns. 


Now when d 4 > 0 this clearly changes the spatio-temporal pattern, as evidenced 
by the sign changes from Case 1 to 2 in Table 1. In this situation the dispersion 
relation is shown in Figure 6 . The resulting patterns are shown Figure 7(a)-(f). 
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Figure 6. Dispersion Plot for two choices of d 4 resembling the two 
cases in Table 1. Notice that when > 0 a band of wavenumbers 
are unstable. 



Remark 6. If r = 0 then (40)-(41) reduces to a model similar to the classical 
predator-prey model with a Holling type II functional response, for which we know 
there cannot occur Turing instability. There is one caveat, the death rate in v now 
becomes nonlinear. So essentially the equations are 


(45) 


d t v = d 2 v xx - I a 2 + w 2 


!)■, 


v + W\- 


uv 


D 1 


(46) d t u = d\u xx + a\u — b 2 u 2 — wo -—. 

u + D 0 

Such systems have been investigated with and without cross and self diffusion [67]. 
In the regular diffusion case, it is known that a Turing instability can exist [56], 
where the nonlinearity in death rate is due to cannibalism. However, to our knowl¬ 
edge, for the specific form of the death rate as above this is not yet known. 

Thus we see that the classical model (4)-(6) can exhibit spatio-temporal patterns, 
apart from just the spatial patterns that were uncovered in [43]. Furthermore addi¬ 
tion of the overcrowding term in (4)-(6), can cause the spatio-temporal patterns to 
change into a purely spatial patterns. It is noted that the inclusion of overcrowding 
with a nonlinear death rate, such as in (45), can lead to Turing instability in the 
classical predator-prey model with Holling type II response. 


6.2. Spatio-temporal chaos. The goal of this section is to investigate spatio- 
temporal chaos in the classical model (4)- ( 6 ). Spatio-temporal chaos is usually 
defined as deterministic dynamics in spatially extended systems that are character¬ 
ized by an apparent randomness in space and time [11]. There is a large literature 









30 


BEAUREGARD, BLACK, PARSHAD, QUANSAH 



i 

0 0.5 1 1.5 2 2.5 3 








Figure 7. The densities of the three species are shown as con¬ 
tour plots in the ad-plane. The long-time simulation yields (a)-(c) 
Turing patterns (d* = 0), that are spatio-temporal, and (d)-(f) 
stripe Turing patterns (dj = 1.6 x 10 _1 ), which are purely spa¬ 
tial. The other parameters are: d\ = 10 -2 , d 2 = 10 -5 , d 3 = 10 -7 , 
ai = 1.79, 02 = 0.8, b 2 = 0.15, c = 0.04, w 3 = 0.55, w\ =2, w 2 = 
0.5, w 3 = 1.2, D 0 = 10, Hi = 13, D 2 = 10, D 3 = 20. 128 grid 
points are used with a temporal step size of .01. 
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on spatio-temporal chaos in PDE, in particular there has been a recent interest on 
spatially extended systems in ecology exhibiting spatio-temporal chaos [31]. How¬ 
ever, most of these works are on two species models, and there is not much literature 
in the three-species case. In [43] we showed diffusion induced temporal chaos in 
(4)-( 6 ), as well as spatial chaos when the domain is enlarged. In [32], various 
patterns non-Turing were uncovered in (4)-(6), in the case of equal diffusion, but 
spatio-temporal chaos was not confirmed. Note, that the appearance of a jagged 
structure in the species density, as seen in [32], which seems to change in time in an 
irregular way, does not necessarily mean that the dynamics is chaotic. One rigoros 
definition of chaos means sensitivity to initial conditions. Thus two initial distribu¬ 
tion, close together, should yield an exponentially growing difference in the species 
distribution at later time. In order to confirm this in (4)-(6), we perform a number 
of tests as in [31]. We run (4)-(6) form a number of different initial conditions, that 
are the same modulo a small perturbation. The parameter set is chosen as in [32]. 
We then look at the difference of the two densities, at each time step in both the 
L°° and L 1 norms. 

Thus we solve (4)-(6) with the following parameter set: d\ = 10~ 5 , d 2 = 
10- 5 , d 3 = 10" 5 , ai = 1.93, a 2 = 1.89, b 2 = 0.06, c = 0.03, w 0 = 1, w x = 
0.5, w 2 = 0.405, u >3 = 1, Do = 10, D\ = 10, D 2 = 10, D 3 = 20. Thus the steady 
state solution for the ODE system is u* = 25, v* = 13, r* = 9. 

The simulations use two different (but close together in L 1 (D),L 2 (fl), L°°(Q) 
norms) initial conditions. The first simulation (which we call r unpe rt) is a pertur¬ 
bation of ( r*,v*,u*) by 0.1cos 2 (x). The second simulation (which we call r pert ) 
is a perturbation of ( r*,v*,u *) by 0.11cos 2 (x). The densities of the species are 
calculated up to the time t = 5000. At each time step in the simulation we compute 

dit') — | \ v unper t (x, t) T per t (x, t) | X 5 

where X = L 1 ( f2), L 2 (fl) and L°°(fl) are used. Then, d(t) is plotted on a log 
scale. In doing son, we observe the exponential growth of the error. This grows at 
an approximate rate of 0.018 > 0. Since this is positive then this is an indicator 
of spatio-temporal chaos. These numerical tests provide experimental evidence to 
the presence of spatio-temporal chaos in the classical model (4)-(6). Figure 8 shows 
the densities of the populations in the xf-plane while Figure 9 gives the error and 
its logarithm till t = 1000 . 


7. Conclusion 

In this work we have proposed and investigated a new model for the control of 
an invasive population. An invasive species population is said to reach catastrophic 
population levels when its population reaches a particular threshold. The mathe¬ 
matical model uses the mathematical construct of finite time blow up which enables 
the model to examine the effect of controls for any particular threshold, especially 
since this level depends on the application. In effect, this construct demonstrates 
that if the invasive population has large enough numbers initially, it can grow to 
explosive levels in finite time: thus wreaking havoc on the ecosystem. Hence, we 
are interested in the influence certain controls will have on the invasive population 
and if its population may be reduced below disastrous levels. This formulation 
yields a clear mathematical problem: Assume that (4)-(6) blows-up in finite time 
for c < w 1 for a given initial condition. Are there controls and features of the 
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Figure 8. The densities of the three species are shown as contour 
plots in the xt -plane for u, v, and r from left to right. The long¬ 
time simulation yields spatio-temporal chaotic patterns. 128 grid 
points are used with a temporal step size of .01. 


model that we can include to modify (4)-(6) so that now there is no blow up in 
the invasive population given the same initial condition? This paper addresses this 
question, suggesting clear controls and improvements to the mathematical model. 
We then investigate these improvements numerically and theoretically. 

In traditional practice, biological control works on the enemy release hypothesis. 
That is releasing an enemy of the invasive species into the ecosystem will lead to 
a decrease in the invasive population. However, non-target effects are prevalent, 
hence this approach is problematic and may create even more devastating impacts 
on the ecosystem [16]. Here, we propose controls that do not use the release of 
biological agents. In particular, we introduce spatial refuges or safe zones for the 
primary food source of the invasive species. Mathematically, this transforms (4)- 
(6) into an indefinite problem [2\ ], that is, where the sign of the coefficient of r 2 
switches between inside and outside of the refuge. We demonstrated numerically 
and in analytically, with some assumptions, that this control can prevent blow up 
and drive the invasive population down. Our numerical experiments suggests that 
there is a delicate balance between the size and location of the refuge and the initial 
condition. In particular, we revealed a logarithmic dependence on the size of the 
initial condition for v versus the critical refuge size to prevent blow-up of r, see 
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Figure 9. The difference d(t) and its logarithm in the L 2 norm for 
the r species, under slightly different initial conditions are shown. 
The difference is seen to grow at an exponential rate of approx¬ 
imately 0.018. Comparable results were found for the L°° and 
L 1 norms. These tests provide experimental evidence that there 
presence of spatio-temporal chaos in the classical model (4)-(6). 


Figure 4. Clearly, the balance is even more pronounced for multiple refuges and in 
higher dimensions. 
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We also improved the mathematical model by incorporating overcrowding ef¬ 
fects, which also may be used a control. We also examined the situation where a 
species may switch its primary food source based in regions of the domain, that 
is, a prey may switch to a predator, hence their roles are reversed. This models 
scenarios where influences on the landscape provide a competitive advantages to 
certain species. Both, role-reversal and overcrowding act as damping mechanisms 
and also may prevent blow-up in the invasive population. In particular, smaller 
refuge sizes in conjunction with role-reversal and overcrowding are required to pre¬ 
vent blow up. Of course, how does one enforce overcrowding in an ecosystem such 
that we obtain this desired effect? Can one devise mechanisms to facilitate this 
dispersal of population? We suggest one approach. Suppose we create a lure, such 
as a pheromone trap, that is placed in the refuge. This would lure the invasive 
species into the patch, where its growth would be controlled. Of course, the species 
would eventually exit the refuge in search for a higher concentration of food. Hence, 
in the future we plan to include spatially dependent diffusion constants to model 
this situation. 

In our mathematical models we confirmed spatio-temporal chaos. We also see 
that overcrowding can effect the sorts of Turing patterns that might form. Since, 
environmental effects are inherently stochastic, part of our future investigations 
introduces stochasticity into the model. It is not known what effect this will have 
on the spatio-temporal chaos or Turing patterns that may emerge. 

It is clear that our new mathematical modeling constructs and results are useful 
analytical and numerical tool for scientists interested in control of invasive species. 
Moreover, the results motivate and encourage numerous avenues of future explo¬ 
ration, many of which are currently under study and will be presented in future 
papers. 
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